rm(list = ls()[!ls()%in%c("drop_attn_fails","n_boot")])
set.seed(221186)
library(tidyverse)
library(readr)
library(psych)
library(survey)
library(cregg)
library(rsample)
library(tictoc)
library(furrr)

# n_boot <- 300

#drop_attn_fails <- TRUE

if(drop_attn_fails){
  
  attn <- "_attn"
  
} else{
  
  attn <- "_full"
  
}

load(paste0("../working/survey_data",attn,".Rdata"))
load(paste0("../working/conjoint_data",attn,".Rdata"))

## Estimate conjoint models

f1 <- mp_selected ~ Descriptive + Substantive + Personalization  + Surrogation + PartisanSurrogation + Responsiveness + Justification  

f1_subrep_fixed <- mp_selected ~ Descriptive  + Personalization  + Surrogation + PartisanSurrogation + Responsiveness + Justification  

f1_partisan_fixed <- mp_selected ~ Descriptive + Substantive + Personalization  + Surrogation + Responsiveness + Justification  

f1_int <- mp_selected ~ 
  Descriptive * Descriptive_FA_Group + 
  Substantive * Substantive_FA_Group + 
  Personalization * Personalization_FA_Group + 
  Surrogation * Surrogation_FA_Group + 
  PartisanSurrogation * Surrogation_FA_Group + 
  Responsiveness * Responsiveness_FA_Group +
  Justification * Justification_FA_Group 

f1_int_linear <- mp_selected ~ 
  Descriptive * desc_fa + 
  Substantive * sub_fa + 
  Personalization * pers_fa + 
  Surrogation * surr_fa + 
  PartisanSurrogation * surr_fa + 
  Responsiveness * resp_fa +
  Justification * just_fa 

## AMCEs

# estimation using cj function for convenience

amces_uk <- cj(uk_long_clean, f1, id = ~ResponseId, weights = uk_long_clean$wgt)
amces_us <- cj(us_long_clean, f1, id = ~ResponseId, weights = us_long_clean$wgt)
amces_de <- cj(de_long_clean, f1, id = ~ResponseId, weights = de_long_clean$wgt)

amces_us <- amces_us %>% 
  mutate(feature = gsub("PartisanSurrogation","Partisan Surrogation",feature),
         feature = gsub("^Surrogation", "Territorial Surrogation", feature),
         feature = factor(feature, levels = c("Descriptive", "Substantive", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")))

amces_uk <- amces_uk %>% 
  mutate(feature = gsub("PartisanSurrogation","Partisan Surrogation",feature),
         feature = gsub("^Surrogation", "Territorial Surrogation", feature),
         feature = factor(feature, levels = c("Descriptive", "Substantive", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")))

amces_de <- amces_de %>% 
  mutate(feature = gsub("PartisanSurrogation","Partisan Surrogation",feature),
         feature = gsub("^Surrogation", "Territorial Surrogation", feature),
         feature = factor(feature, levels = c("Descriptive", "Substantive", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")))

xlims <- range(c(amces_uk$lower,amces_uk$upper, amces_us$lower, amces_us$upper, amces_de$lower, amces_de$upper), na.rm = TRUE)

pdf(paste0("../output/plots/conjoint_amce/conjoint_amce_uk",attn,".pdf"),6,7)
amce_plot_uk <- ggplot(amces_uk, aes(x = estimate, xmax = upper, xmin = lower, y = level, group = feature)) + 
  geom_point() + geom_pointrange() + 
  facet_wrap(~feature, scale = "free_y", ncol = 1) + 
  theme_bw() + geom_vline(xintercept = 0, linetype = 2) + 
  ylab("") + xlab("AMCE") + 
  ggtitle("UK") + 
  xlim(xlims)
print(amce_plot_uk)
dev.off()

pdf(paste0("../output/plots/conjoint_amce/conjoint_amce_us",attn,".pdf"),6,7)
amce_plot_us <- ggplot(amces_us, aes(x = estimate, xmax = upper, xmin = lower, y = level, group = feature)) + 
  geom_point() + geom_pointrange() + 
  facet_wrap(~feature, scale = "free_y", ncol = 1) + 
  theme_bw() + geom_vline(xintercept = 0, linetype = 2) + 
  ylab("") + xlab("AMCE") + 
  ggtitle("US") + 
  xlim(xlims)
print(amce_plot_us)
dev.off()

pdf(paste0("../output/plots/conjoint_amce/conjoint_amce_de",attn,".pdf"),6,7)
amce_plot_de <- ggplot(amces_de, aes(x = estimate, xmax = upper, xmin = lower, y = level, group = feature)) + 
  geom_point() + geom_pointrange() + 
  facet_wrap(~feature, scale = "free_y", ncol = 1) + 
  theme_bw() + geom_vline(xintercept = 0, linetype = 2) + 
  ylab("") + xlab("AMCE") + 
  ggtitle("Germany") + 
  xlim(xlims)
print(amce_plot_de)
dev.off()

amces_de$country <- "Germany"
amces_uk$country <- "UK"
amces_us$country <- "US"

amces_all <- rbind(amces_de, amces_uk, amces_us)

levels(amces_all$level) <- gsub("district|constituency","district", levels(amces_all$level))
levels(amces_all$level) <- gsub("MP","Representative", levels(amces_all$level))
levels(amces_all$level) <- gsub("Representative","Rep.", levels(amces_all$level))
amces_all$feature <- gsub("Territorial ", "Territorial \n", amces_all$feature)
amces_all$feature <- gsub("Partisan ", "Partisan \n", amces_all$feature)

amces_all$feature <- factor(amces_all$feature, levels = c("Substantive", "Descriptive", "Territorial \nSurrogation", "Partisan \nSurrogation", "Justification", "Personalization", "Responsiveness"))

ggplot(amces_all, aes(x = estimate, xmax = upper, xmin = lower, y = level, group = feature)) + 
  geom_point() + geom_pointrange() + 
  facet_grid(feature ~ country, scales = "free_y" ) + 
  theme_bw() + geom_vline(xintercept = 0, linetype = 2) + 
  ylab("") + xlab("AMCE") + 
  theme(strip.background = element_blank(), strip.placement = "outside")

ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined",attn,".pdf"), width = 16, height = 7.5)

# Multiple corrections adjustment

amces_all %>%
  mutate(p.new = p.adjust(p, "BH"),
         sig = p.new <= 0.05,
         sig = ifelse(is.na(sig),FALSE, sig),
         sig = ifelse(sig, "<= 0.05", "> 0.05")) %>%
  ggplot(aes(x = estimate, xmax = upper, xmin = lower, y = level, group = feature, col = sig)) + 
  geom_point() + geom_pointrange() + 
  facet_grid(feature ~ country, scales = "free_y" ) + 
  theme_bw() + geom_vline(xintercept = 0, linetype = 2) + 
  ylab("") + xlab("AMCE") + 
  theme(strip.background = element_blank(), strip.placement = "outside") +     
  scale_color_manual("Adjusted p-value",values = c("black","gray")) +
  theme(legend.position = "bottom") + 
  guides(colour = guide_legend(title.position="top", title.hjust = 0.5)) 
ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined",attn,"_BH.pdf"), width = 16, height = 7.5)

x <- amces_all %>%
  mutate(p.new = p.adjust(p, "BH"))
x[x$feature == "Descriptive" & x$country=="Germany",]
x[which(x$p <= 0.05 & (! x$p.new <= 0.05)),]

## Create table of outputs for pedants at APSR

sink(paste0("../output/tables/amces",attn,".tex"))

if(attn == "_full") {
  title <- "Average marginal component effects (Full sample)"
} else{
  title <- "Average marginal component effects (Weighted sample)"  
}


amces_all %>%
  mutate(est = paste0(round(estimate,2), " (", round(std.error,2), ")"),
         est = ifelse(grepl("NA", est),"-", est)) %>%
  select(feature, level, est, country) %>%
  pivot_wider(names_from = country,
              values_from = est) %>%
  mutate(feature = as.character(feature),
         feature = gsub("\\n", "", feature),
         feature = factor(feature, levels = c("Substantive", "Descriptive", "Personalization", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Responsiveness"))) %>%
  arrange(feature) %>%
  mutate(level = trimws(as.character(level))) %>%
  rename(Level = level,
         Variable = feature) %>%
  kableExtra::kable(format = "latex", booktabs = TRUE, escape = TRUE, align = "lllllll", linesep = "", caption = title, label = paste0("amces_table",attn)) %>%
  kableExtra::kable_styling(font_size = 10) %>%
  kableExtra::add_footnote(paste0("Sample sizes: Germany = ", nrow(de_long_clean),", UK = ",nrow(uk_long_clean),", US = ", nrow(us_long_clean),". The dependent variable is equal to 1 if an MP was selected as the better representative in a given comparison and 0 otherwise. Respondent-clustered standard errors in parentheses."), notation = "number")
sink()

## ############
## Do conjoint estimates vary by task number?
## ############

get_amces_by_task <- function(task_number){
  
  # estimation using cj function for convenience
  
  uk_tmp <- uk_long_clean[uk_long_clean$task == task_number,]
  us_tmp <- us_long_clean[us_long_clean$task == task_number,]
  de_tmp <- de_long_clean[de_long_clean$task == task_number,]
  
  amces_uk <- cj(uk_tmp, f1, id = ~ResponseId, weights = uk_tmp$wgt)
  amces_us <- cj(us_tmp, f1, id = ~ResponseId, weights = us_tmp$wgt)
  amces_de <- cj(de_tmp, f1, id = ~ResponseId, weights = de_tmp$wgt)
  
  amces_us <- amces_us %>% 
    mutate(feature = gsub("PartisanSurrogation","Partisan Surrogation",feature),
           feature = gsub("^Surrogation", "Territorial Surrogation", feature),
           feature = factor(feature, levels = c("Descriptive", "Substantive", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")))
  
  amces_uk <- amces_uk %>% 
    mutate(feature = gsub("PartisanSurrogation","Partisan Surrogation",feature),
           feature = gsub("^Surrogation", "Territorial Surrogation", feature),
           feature = factor(feature, levels = c("Descriptive", "Substantive", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")))
  
  amces_de <- amces_de %>% 
    mutate(feature = gsub("PartisanSurrogation","Partisan Surrogation",feature),
           feature = gsub("^Surrogation", "Territorial Surrogation", feature),
           feature = factor(feature, levels = c("Descriptive", "Substantive", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")))
  
  amces_de$country <- "Germany"
  amces_uk$country <- "UK"
  amces_us$country <- "US"
  
  amces_all <- rbind(amces_de, amces_uk, amces_us)
  
  levels(amces_all$level) <- gsub("district|constituency","district", levels(amces_all$level))
  levels(amces_all$level) <- gsub("MP","Representative", levels(amces_all$level))
  levels(amces_all$level) <- gsub("Representative","Rep.", levels(amces_all$level))
  amces_all$feature <- gsub("Territorial ", "Territorial \n", amces_all$feature)
  amces_all$feature <- gsub("Partisan ", "Partisan \n", amces_all$feature)
  
  amces_all$feature <- factor(amces_all$feature, levels = c("Substantive", "Descriptive", "Territorial \nSurrogation", "Partisan \nSurrogation", "Justification", "Personalization", "Responsiveness"))
  
  amces_all$task <- task_number
  return(amces_all)
  
}

amces_by_task <- lapply(1:5, function(i) get_amces_by_task(i))

bind_rows(amces_by_task) %>%
  ggplot(aes(x = estimate, xmax = upper, xmin = lower, y = level, group = feature, col = as.factor(task))) + 
  geom_point(position = position_dodge2(.4)) + geom_pointrange(position = position_dodge2(.4)) + 
  facet_grid(feature ~ country, scales = "free_y" ) + 
  theme_bw() + geom_vline(xintercept = 0, linetype = 2) + 
  ylab("") + xlab("AMCE") + 
  theme(strip.background = element_blank(), strip.placement = "outside") + 
  labs(color = "Conjoint task number") + 
  theme(legend.position = "bottom")

ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined_by_task",attn,".pdf"), width = 16, height = 7.5)

get_tables_by_task <- function(this_country, country_data){

  tmp <- bind_rows(amces_by_task) %>%
    filter(country == this_country) %>%
    mutate(est = paste0(round(estimate,3), " (", round(std.error,3),")"),
           est = ifelse(grepl("NA", est), "-", est)) %>%
    select(feature, level, est, task) %>%
    mutate(level = trimws(as.character(level)),
           feature = gsub("\\n", "", feature)) %>%
    pivot_wider(values_from = est, names_from = task)
    
  names(tmp) <- c("Dimension", "Level", paste0("Task ", 1:5))
  tmp %>%
    kableExtra::kable(format = "latex", booktabs = TRUE, escape = TRUE, align = "lllllll", linesep = "", caption = paste0(paste0("Average marginal component effects, ", this_country, ", by task number.")),
                      label = paste0("amce_by_task_",this_country)) %>%
    kableExtra::kable_styling(font_size = 10) %>%
    kableExtra::add_footnote(paste0("Sample size: ", nrow(country_data),". The dependent variable is equal to 1 if an MP was selected as the better representative in a given comparison and 0 otherwise. Respondent-clustered standard errors in parentheses."), notation = "number") %>%
    print()
}

sink("../output/tables/conjoint/by_task_de.tex")
get_tables_by_task("Germany", de_long_clean)
sink()

sink("../output/tables/conjoint/by_task_us.tex")
get_tables_by_task("US", us_long_clean)
sink()

sink("../output/tables/conjoint/by_task_uk.tex")
get_tables_by_task("UK", uk_long_clean)
sink()

## ############
## Conditional AMCEs by Substantive Representation, Party Representation, Gender, Race and Sexuality
## ############

arbitrary_subgroup_amces <- function(this_model = f1, variable_name = "Partisan Surrogation"){
  
  these_groups <- levels(uk_long_clean$tmp_group)
  
  uk_long_clean_g1 <- uk_long_clean %>% filter(tmp_group == these_groups[1])
  uk_long_clean_g2 <- uk_long_clean %>% filter(tmp_group == these_groups[2])
  
  us_long_clean_g1 <- us_long_clean %>% filter(tmp_group == these_groups[1])
  us_long_clean_g2 <- us_long_clean %>% filter(tmp_group == these_groups[2])
  
  de_long_clean_g1 <- de_long_clean %>% filter(tmp_group == these_groups[1])
  de_long_clean_g2 <- de_long_clean %>% filter(tmp_group == these_groups[2])
  
  all_clean <- rbind(uk_long_clean[,c("ResponseId", "tmp_group", "wgt", all.vars(f1))],
        us_long_clean[,c("ResponseId","tmp_group", "wgt", all.vars(f1))],
        de_long_clean[,c("ResponseId", "tmp_group", "wgt", all.vars(f1))])
  
  levels(all_clean$Surrogation) <- gsub("district|constituency","district", levels(all_clean$Surrogation))
  levels(all_clean$Surrogation) <- gsub("MP","Representative", levels(all_clean$Surrogation))
  levels(all_clean$Surrogation) <- gsub("Representative","Rep.", levels(all_clean$Surrogation))
  
  levels(all_clean$PartisanSurrogation) <- gsub("district|constituency","district", levels(all_clean$PartisanSurrogation))
  levels(all_clean$PartisanSurrogation) <- gsub("MP","Representative", levels(all_clean$PartisanSurrogation))
  levels(all_clean$PartisanSurrogation) <- gsub("Representative","Rep.", levels(all_clean$PartisanSurrogation))
  
  all_clean_g1 <- all_clean %>% filter(tmp_group == these_groups[1])
  all_clean_g2 <- all_clean %>% filter(tmp_group == these_groups[2])
  
  amces_uk_g1 <- cj(uk_long_clean_g1, this_model, 
                    id = ~ResponseId, weights = uk_long_clean_g1$wgt)
  amces_uk_g2 <- cj(uk_long_clean_g2, this_model, 
                    id = ~ResponseId, weights = uk_long_clean_g2$wgt)
  
  amces_us_g1 <- cj(us_long_clean_g1, this_model, 
                    id = ~ResponseId, weights = us_long_clean_g1$wgt)
  amces_us_g2 <- cj(us_long_clean_g2, this_model, 
                    id = ~ResponseId, weights = us_long_clean_g2$wgt)
  
  amces_de_g1 <- cj(de_long_clean_g1, this_model, 
                    id = ~ResponseId, weights = de_long_clean_g1$wgt)
  amces_de_g2 <- cj(de_long_clean_g2, this_model, 
                    id = ~ResponseId, weights = de_long_clean_g2$wgt)
  
  amces_all_g1 <- cj(all_clean_g1, this_model, 
                    id = ~ResponseId, weights = all_clean_g1$wgt)
  amces_all_g2 <- cj(all_clean_g2, this_model, 
                    id = ~ResponseId, weights = all_clean_g2$wgt)
  
  amces_all_g1$group <- amces_uk_g1$group <- amces_us_g1$group <- amces_de_g1$group <- these_groups[1]
  amces_all_g2$group <- amces_uk_g2$group <- amces_us_g2$group <- amces_de_g2$group <- these_groups[2]
  
  amces_uk_g2$country <- amces_uk_g1$country <- "UK"
  amces_us_g2$country <- amces_us_g1$country <- "US"
  amces_de_g2$country <- amces_de_g1$country <- "DE"
  amces_all_g2$country <- amces_all_g1$country <- "Pooled"
  
  amces_all <- rbind(amces_uk_g2, amces_uk_g1,
                     amces_us_g2, amces_us_g1,
                     amces_de_g2, amces_de_g1,
                     amces_all_g2, amces_all_g1)
  
  amces_all <- tibble(amces_all) %>% 
    mutate(feature = gsub("PartisanSurrogation", "Partisan Surrogation", feature),
           feature = gsub("^Surrogation", "Territorial Surrogation", feature),
           feature = factor(feature, levels = c("Substantive", "Descriptive", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")),
           country = factor(country, levels = c("Pooled", "UK", "US", "DE"))) 
  
  levels(amces_all$level) <- gsub("district|constituency","district", levels(amces_all$level))
  levels(amces_all$level) <- gsub("MP","Representative", levels(amces_all$level))
  levels(amces_all$level) <- gsub("Representative","Rep.", levels(amces_all$level))
  levels(amces_all$feature) <- gsub("Territorial","Territorial \n", levels(amces_all$feature))
  levels(amces_all$feature) <- gsub("Partisan","Partisan \n", levels(amces_all$feature))
  
  amces_all$upper[is.na(amces_all$upper)] <- 0
  amces_all$lower[is.na(amces_all$lower)] <- 0
  
  amces_tabs <- amces_all %>%
    mutate(est = paste0(round(estimate, 3), " (", round(std.error,3),")"),
           est = ifelse(grepl("NA", est), "-", est)) %>%
    filter(est != "-") %>%
    select(feature, level, group, est, country) %>%
    pivot_wider(values_from = est,
                names_from = country) 
  
  amces_tabs <- amces_tabs %>% arrange(group, feature, level) %>%
    mutate(feature = gsub("\\n", "", as.character(feature)),
           level = trimws(as.character(level))) %>%
    rename(Dimension = feature,
           Level = level)
  
  amces_tabs <- amces_tabs %>%
    select(-group) %>%
    kableExtra::kable(format = "latex", booktabs = TRUE, escape = TRUE, 
                      align = "lllllll", linesep = "", 
                      caption = paste0("Average Marginal Component Effects by ", variable_name, "."), 
                      label = paste0("amce_by_",variable_name)) %>%
    kableExtra::kable_styling() %>%
    kableExtra::pack_rows(index = table(fct_inorder(amces_tabs$group)))%>%
    kableExtra::kable_styling(font_size = 10) %>%
    kableExtra::add_footnote(paste0("Sample sizes: Germany = ", nrow(de_long_clean),", UK = ",nrow(uk_long_clean),", US = ",nrow(us_long_clean),". The dependent variable is equal to 1 if an MP was selected as the better representative in a given comparison and 0 otherwise. Respondent-clustered standard errors in parentheses."), notation = "number")
  
  p1 <- amces_all %>%
    ggplot(aes(x = estimate, xmax = upper, xmin = lower, y = level, group = group, col = group)) +
    geom_pointrange(position = position_dodge(.5)) + 
    facet_grid(feature ~ country, scales = "free_y" ) + 
    theme_bw() + geom_vline(xintercept = 0, linetype = 2) + 
    ylab("") + xlab("AMCE") + 
    theme(strip.background = element_blank(), strip.placement = "outside",
          legend.position = "bottom") + 
    scale_color_manual("", values = c("black", "gray")) 
  
  return(list(p1 = p1, amces_tabs = amces_tabs))
  
}



## AMCEs by Partisan Surrogation

f1_no_partisan <- mp_selected ~ Descriptive + Substantive  + Personalization  + Surrogation  + Responsiveness + Justification  

uk_long_clean <- uk_long_clean %>% 
  mutate(tmp_group = PartisanSurrogation == "MP from respondent's party",
         tmp_group = factor(ifelse(tmp_group, "Respondent's party", "Another Party")))

us_long_clean <- us_long_clean %>% 
  mutate(tmp_group = PartisanSurrogation == "Representative from respondent's party",
         tmp_group = factor(ifelse(tmp_group, "Respondent's party", "Another Party")))

de_long_clean <- de_long_clean %>% 
  mutate(tmp_group = PartisanSurrogation == "Representative from respondent's party",
         tmp_group = factor(ifelse(tmp_group, "Respondent's party", "Another Party")))

partrep_amces <- arbitrary_subgroup_amces(this_model = f1_no_partisan, variable_name = "Partisan Surrogation")

print(partrep_amces$p1)
ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined_partrep_conditional",attn,".pdf"), width = 12, height = 8)

sink(paste0("../output/tables/conjoint/conjoint_amce_combined_partrep_conditional",attn,".tex"))
print(partrep_amces$amces_tabs)
sink()

## AMCEs by Substantive Representation

f1_no_substantive <- mp_selected ~ Descriptive  + Personalization  + Surrogation + PartisanSurrogation + Responsiveness + Justification  

uk_long_clean <- uk_long_clean %>% mutate(tmp_group = Substantive == "    Congruent",
                                          tmp_group = factor(ifelse(tmp_group, "Congruent", "Incongruent")))
us_long_clean <- us_long_clean %>% mutate(tmp_group = Substantive == "    Congruent",
                                          tmp_group = factor(ifelse(tmp_group, "Congruent", "Incongruent")))
de_long_clean <- de_long_clean %>% 
  mutate(tmp_group = Substantive == "    Congruent",
         tmp_group = factor(ifelse(tmp_group, "Congruent", "Incongruent")))


substantiverep_amces <- arbitrary_subgroup_amces(this_model = f1_no_substantive, variable_name = "Substantive Representation")

print(substantiverep_amces$p1)

ggsave("../output/plots/conjoint_amce/conjoint_amce_combined_subrep_conditional_attn.pdf", width = 12, height = 8)

sink(paste0("../output/tables/conjoint/conjoint_amce_combined_subrep_conditional",attn,".tex"))
print(substantiverep_amces$amces_tabs)
sink()

## AMCEs by Gender

uk_long_clean <- uk_long_clean %>% mutate(tmp_group = Gender)
us_long_clean <- us_long_clean %>% mutate(tmp_group = Gender)
de_long_clean <- de_long_clean %>% mutate(tmp_group = factor(ifelse(Gender == "Männlich", "Male", "Female"), levels = levels(uk_long_clean$tmp_group)))

gender_amces <- arbitrary_subgroup_amces(this_model = f1, variable_name = "Gender")

print(gender_amces$p1)
ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined_gender_conditional",attn,".pdf"), width = 12, height = 8)

sink(paste0("../output/tables/conjoint/conjoint_amce_combined_gender_conditional",attn,".tex"))
print(gender_amces$amces_tabs)
sink()



## AMCEs by Ethnicity

uk_long_clean <- uk_long_clean %>% mutate(tmp_group = factor(ifelse(ethnicity == "White", "White/No Migrant Background", "Other")))
us_long_clean <- us_long_clean %>% mutate(tmp_group = factor(ifelse(ethnicity == "White", "White/No Migrant Background", "Other")))
de_long_clean <- de_long_clean %>% mutate(tmp_group = factor(ifelse(ethnicity == "No migrant background", "White/No Migrant Background", "Other")))

ethnicity_amces <- arbitrary_subgroup_amces(this_model = f1, variable_name = "Race/Ethnicity/Migratory Background")

print(ethnicity_amces$p1)
ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined_race_conditional",attn,".pdf"), width = 12, height = 8)

sink(paste0("../output/tables/conjoint/conjoint_amce_combined_ethnicity_conditional",attn,".tex"))
print(ethnicity_amces$amces_tabs)
sink()

## AMCEs by Sexuality

uk_long_clean <- uk_long_clean %>% mutate(tmp_group = factor(ifelse(Sexuality == "Heterosexual", "Heterosexual", "Other")))
us_long_clean <- us_long_clean %>% mutate(tmp_group = factor(ifelse(Sexuality == "Heterosexual", "Heterosexual", "Other")))
de_long_clean <- de_long_clean %>% mutate(tmp_group = factor(ifelse(Sexuality == "Heterosexuell", "Heterosexual", "Other")))

sex_amces <- arbitrary_subgroup_amces(this_model = f1, variable_name = "Sexuality")

print(sex_amces$p1)
ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined_sexuality_conditional",attn,".pdf"), width = 12, height = 8)

sink(paste0("../output/tables/conjoint/conjoint_amce_combined_sex_conditional",attn,".tex"))
print(sex_amces$amces_tabs)
sink()



## AMCEs by seeing conjoint first/seeing conjoint second

uk_long_clean <- uk_long_clean %>% mutate(tmp_group = factor(ifelse(saw_conjoint_first, "Conjoint First", "Batteries First")))
us_long_clean <- us_long_clean %>% mutate(tmp_group = factor(ifelse(saw_conjoint_first, "Conjoint First", "Batteries First")))
de_long_clean <- de_long_clean %>% mutate(tmp_group = factor(ifelse(saw_conjoint_first, "Conjoint First", "Batteries First")))

amce_randomization_order <- arbitrary_subgroup_amces(this_model = f1)
amce_randomization_order$p1
ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined_randomization_order_conditional",attn,".pdf"), width = 12, height = 8)

## Decomposing descriptive representation effects 

uk_long_clean$Descriptive_original <- uk_long_clean$Descriptive
us_long_clean$Descriptive_original <- us_long_clean$Descriptive
de_long_clean$Descriptive_original <- de_long_clean$Descriptive

uk_long_clean <- uk_long_clean %>% mutate(Descriptive = factor((gender_match == descriptive_1_mp | gender_match == descriptive_2_mp), labels = c("Different gender", "Same Gender")))
us_long_clean <- us_long_clean %>% mutate(Descriptive = factor((gender_match == descriptive_1_mp | gender_match == descriptive_2_mp), labels = c("Different gender", "Same Gender")))
de_long_clean <- de_long_clean %>% mutate(Descriptive = factor((gender_match == descriptive_1_mp | gender_match == descriptive_2_mp), labels = c("Different gender", "Same Gender")))

uk_long_clean <- uk_long_clean %>% mutate(tmp_group = Gender)
us_long_clean <- us_long_clean %>% mutate(tmp_group = Gender)
de_long_clean <- de_long_clean %>% mutate(tmp_group = factor(ifelse(Gender == "Männlich", "Male", "Female"), levels = levels(uk_long_clean$tmp_group)))

amce_separate_descriptive <- arbitrary_subgroup_amces(this_model = f1)

print(amce_separate_descriptive$p1)
ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined_gender_conditional_separate_descriptive_effects",attn,".pdf"), width = 16, height = 7.5)

uk_long_clean <- uk_long_clean %>% mutate(Descriptive = factor((ethnicity_match == descriptive_1_mp | ethnicity_match == descriptive_2_mp), labels = c("Different ethnicity/race", "Same ethnicity/race")))
us_long_clean <- us_long_clean %>% mutate(Descriptive = factor((ethnicity_match == descriptive_1_mp | ethnicity_match == descriptive_2_mp), labels = c("Different ethnicity/race", "Same ethnicity/race")))
de_long_clean <- de_long_clean %>% mutate(Descriptive = factor((ethnicity_match == descriptive_1_mp | ethnicity_match == descriptive_2_mp), labels = c("Different ethnicity/race", "Same ethnicity/race")))

uk_long_clean <- uk_long_clean %>% mutate(tmp_group = factor(ifelse(grepl("White", ethnicity), "White/no migrant background", "Other race/ethnicity")))
us_long_clean <- us_long_clean %>% mutate(tmp_group = factor(ifelse(grepl("White", ethnicity), "White/no migrant background", "Other race/ethnicity")))
de_long_clean <- de_long_clean %>% mutate(tmp_group = factor(ifelse(grepl("Nein", Migr_back), "White/no migrant background", "Other race/ethnicity")))

amce_separate_descriptive <- arbitrary_subgroup_amces(this_model = f1)

print(amce_separate_descriptive$p1)

ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined_race_conditional_separate_descriptive_effects",attn,".pdf"), width = 16, height = 7.5)


uk_long_clean <- uk_long_clean %>% mutate(Descriptive = factor((sexuality_match == descriptive_1_mp | sexuality_match == descriptive_2_mp), labels = c("Different sexuality", "Same sexuality")))
us_long_clean <- us_long_clean %>% mutate(Descriptive = factor((sexuality_match == descriptive_1_mp | sexuality_match == descriptive_2_mp), labels = c("Different sexuality", "Same sexuality")))
de_long_clean <- de_long_clean %>% mutate(Descriptive = factor((sexuality_match == descriptive_1_mp | sexuality_match == descriptive_2_mp), labels = c("Different sexuality", "Same sexuality")))

uk_long_clean <- uk_long_clean %>% mutate(tmp_group = factor(ifelse(grepl("heterosexual", sexuality_match), "Heterosexual", "Other sexuality")))
us_long_clean <- us_long_clean %>% mutate(tmp_group = factor(ifelse(grepl("heterosexual", sexuality_match), "Heterosexual", "Other sexuality")))
de_long_clean <- de_long_clean %>% mutate(tmp_group = factor(ifelse(grepl("heterosexuell", sexuality_match), "Heterosexual", "Other sexuality")))

amce_separate_descriptive <- arbitrary_subgroup_amces(this_model = f1)

print(amce_separate_descriptive$p1)

ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined_sexuality_conditional_separate_descriptive_effects",attn,".pdf"), width = 16, height = 7.5)

uk_long_clean$Descriptive <- uk_long_clean$Descriptive_original
us_long_clean$Descriptive <- us_long_clean$Descriptive_original
de_long_clean$Descriptive <- de_long_clean$Descriptive_original


## Linear and binned groups of factor analysis measures

svydesign_out_uk <- survey::svydesign(ids = ~ResponseId, weights = uk_long_clean$wgt, data = uk_long_clean)
uk_conjoint_interaction <- survey::svyglm(f1_int, design = svydesign_out_uk)
uk_conjoint_interaction_linear <- survey::svyglm(f1_int_linear, design = svydesign_out_uk)

# US

svydesign_out_us <- survey::svydesign(ids = ~ResponseId, weights = us_long_clean$wgt, data = us_long_clean)
us_conjoint_interaction <- survey::svyglm(f1_int, design = svydesign_out_us)
us_conjoint_interaction_linear <- survey::svyglm(f1_int_linear, design = svydesign_out_us)

# DE

svydesign_out_de <- survey::svydesign(ids = ~ResponseId, weights = de_long_clean$wgt, data = de_long_clean)
de_conjoint_interaction <- survey::svyglm(f1_int, design = svydesign_out_de)
de_conjoint_interaction_linear <- survey::svyglm(f1_int_linear, design = svydesign_out_de)

## Are the interaction effects positive and significant for each dimension?

uk_coefs <- broom::tidy(uk_conjoint_interaction_linear, conf.int = TRUE)
us_coefs <- broom::tidy(us_conjoint_interaction_linear, conf.int = TRUE)
de_coefs <- broom::tidy(de_conjoint_interaction_linear, conf.int = TRUE)

uk_coefs_group <- broom::tidy(uk_conjoint_interaction, conf.int = TRUE)
us_coefs_group <- broom::tidy(us_conjoint_interaction, conf.int = TRUE)
de_coefs_group <- broom::tidy(de_conjoint_interaction, conf.int = TRUE)

## Bootstrap intervals

acme_interaction_boot <- function(x, my_model = f1_int_linear, base_data = uk_long_clean){
  
  x <- x %>%
    select(c(ResponseId, fa_names)) %>%
    right_join(base_data[,!names(base_data) %in% fa_names], by = "ResponseId", relationship = "many-to-many")
  x <- survey::svydesign(ids = ~ResponseId, weights = x$wgt, data = x)
  mod_out <- survey::svyglm(my_model, design = x)  
  coef(mod_out)
  
}

source("est_factors_v2.R")

out_boot <- est_factors(n_boots = n_boot, attn = attn)

fa_names <- c("pers_fa", "sub_fa", "just_fa", "resp_fa", "desc_fa", "surr_fa")

# Linear specification for factors

boot_uk <- out_boot$boots_uk %>%
  mutate(coef = map(out, acme_interaction_boot, f1_int_linear, uk_long_clean)) %>%
  select(coef) %>% unlist() %>%
  matrix(ncol = n_boot, byrow = FALSE)

boot_us <- out_boot$boots_us %>%
  mutate(coef = map(out, acme_interaction_boot, f1_int_linear, us_long_clean)) %>%
  select(coef) %>% unlist() %>%
  matrix(ncol = n_boot, byrow = FALSE)

boot_de <- out_boot$boots_de %>%
  mutate(coef = map(out, acme_interaction_boot, f1_int_linear, de_long_clean)) %>%
  select(coef) %>% unlist() %>%
  matrix(ncol = n_boot, byrow = FALSE)

rownames(boot_uk) <- names(coef(uk_conjoint_interaction_linear))
rownames(boot_us) <- names(coef(us_conjoint_interaction_linear))
rownames(boot_de) <- names(coef(de_conjoint_interaction_linear))

# Binned specification for factors

boot_uk_group <- out_boot$boots_uk %>%
  mutate(coef = map(out, acme_interaction_boot, f1_int, uk_long_clean)) %>%
  select(coef) %>% unlist() %>%
  matrix(ncol = n_boot, byrow = FALSE)

boot_us_group <- out_boot$boots_us %>%
  mutate(coef = map(out, acme_interaction_boot, f1_int, us_long_clean)) %>%
  select(coef) %>% unlist() %>%
  matrix(ncol = n_boot, byrow = FALSE)

boot_de_group <- out_boot$boots_de %>%
  mutate(coef = map(out, acme_interaction_boot, f1_int, de_long_clean)) %>%
  select(coef) %>% unlist() %>%
  matrix(ncol = n_boot, byrow = FALSE)

rownames(boot_uk_group) <- names(coef(uk_conjoint_interaction))
rownames(boot_us_group) <- names(coef(us_conjoint_interaction))
rownames(boot_de_group) <- names(coef(de_conjoint_interaction))

## Combine with point estimates

uk_coefs$boot_high <- apply(boot_uk,1,quantile, 0.975)
us_coefs$boot_high <- apply(boot_us,1,quantile, 0.975)
de_coefs$boot_high <- apply(boot_de,1,quantile, 0.975)

uk_coefs$boot_low <- apply(boot_uk,1,quantile, 0.025)
us_coefs$boot_low <- apply(boot_us,1,quantile, 0.025)
de_coefs$boot_low <- apply(boot_de,1,quantile, 0.025)

uk_coefs_group$boot_high <- apply(boot_uk_group,1,quantile, 0.975)
us_coefs_group$boot_high <- apply(boot_us_group,1,quantile, 0.975)
de_coefs_group$boot_high <- apply(boot_de_group,1,quantile, 0.975)

uk_coefs_group$boot_low <- apply(boot_uk_group,1,quantile, 0.025)
us_coefs_group$boot_low <- apply(boot_us_group,1,quantile, 0.025)
de_coefs_group$boot_low <- apply(boot_de_group,1,quantile, 0.025)


# Plot interaction effects

## Linear effects

uk_interaction_ests <- uk_coefs %>% filter(grepl(":",term)) %>%
  mutate(dimension = case_when(grepl("Descriptive1", term) ~ "Descriptive (Partial)",
                               grepl("Descriptive2", term) ~ "Descriptive (Full)",
                               grepl("Substantive", term) ~ "Substantive",
                               grepl("Personalization", term) ~ "Personalization",
                               grepl("SurrogationMP for another constituency", term) ~ "Territorial Surrogation",
                               grepl("PartisanSurrogationMP", term) ~ "Partisan Surrogation",
                               grepl("Responsiveness", term) ~ "Responsiveness",
                               grepl("Justification", term) ~ "Justification")) 

us_interaction_ests <- us_coefs %>% filter(grepl(":",term)) %>%
  mutate(dimension = case_when(grepl("Descriptive1", term) ~ "Descriptive (Partial)",
                               grepl("Descriptive2", term) ~ "Descriptive (Full)",
                               grepl("Substantive", term) ~ "Substantive",
                               grepl("Personalization", term) ~ "Personalization",
                               grepl("SurrogationRepresentative for another district", term) ~ "Territorial Surrogation",
                               grepl("PartisanSurrogation", term) ~ "Partisan Surrogation",
                               grepl("Responsiveness", term) ~ "Responsiveness",
                               grepl("Justification", term) ~ "Justification"))

de_interaction_ests <- de_coefs %>% filter(grepl(":",term)) %>%
  mutate(dimension = case_when(grepl("Descriptive1", term) ~ "Descriptive (Partial)",
                               grepl("Descriptive2", term) ~ "Descriptive (Full)",
                               grepl("Substantive", term) ~ "Substantive",
                               grepl("Personalization", term) ~ "Personalization",
                               grepl("SurrogationMP for another constituency in another", term) ~ "Territorial Surrogation", # Note this is only one of the DE Territorial surrogation dimensions
                               grepl("PartisanSurrogation", term) ~ "Partisan Surrogation",
                               grepl("Responsiveness", term) ~ "Responsiveness",
                               grepl("Justification", term) ~ "Justification")) %>%
  filter(!is.na(dimension))

uk_interaction_ests$country <- "UK"
us_interaction_ests$country <- "US"
de_interaction_ests$country <- "Germany"

interaction_linear_out <- rbind(uk_interaction_ests, us_interaction_ests, de_interaction_ests) %>%
  mutate(dimension = factor(dimension, levels = rev(c("Substantive", "Descriptive (Partial)", "Descriptive (Full)", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")))) 

interaction_linear_out %>%
  ggplot(aes(x = estimate, y = dimension, group = country)) + 
  geom_pointrange(aes(xmin = boot_low, xmax = boot_high)) + 
  theme_bw() + 
  geom_vline(xintercept = 0, linetype = 2) + 
  facet_grid(~country) + 
  xlab("Average Marginal Component Interaction Effect") + 
  ylab("")

ggsave(paste0("../output/plots/conjoint_amce_validation/amce_interaction_linear",attn,".pdf"), width = 12, height = 6)

sink(paste0("../output/tables/conjoint/amce_interaction_linear",attn,".tex"))
interaction_linear_out %>%
  mutate(est = paste0(round(estimate, 3), " (",round(std.error,3), ")")) %>%
  select(dimension, est, country) %>%
  pivot_wider(names_from = country,
              values_from = est) %>%
  arrange(desc(dimension)) %>%
  kableExtra::kable(format = "latex", booktabs = TRUE, escape = TRUE, 
                    align = "lllllll", linesep = "", 
                    caption = paste0("Average Marginal Component Interaction Effects (Linear)."), 
                    label = paste0("amce_interaction_effects_linear",attn)) %>%
  kableExtra::kable_styling() %>%
  kableExtra::kable_styling(font_size = 10) %>%
  kableExtra::add_footnote(paste0("Sample sizes: Germany = ", nrow(de_long_clean),", UK = ",nrow(uk_long_clean),", US = ", nrow(us_long_clean),". The dependent variable is equal to 1 if an MP was selected as the better representative in a given comparison and 0 otherwise. Respondent-clustered standard errors in parentheses."), notation = "number") %>%
  print()
sink()

## Group effects

uk_interaction_ests_group <- uk_coefs_group %>% filter(grepl(":",term)) %>%
  mutate(group = case_when(grepl("GroupMid", term) ~ "Mid",
                           grepl("GroupHi", term) ~ "High"),
         dimension = case_when(grepl("Descriptive1", term) ~ "Descriptive (Partial)",
                               grepl("Descriptive2", term) ~ "Descriptive (Full)",
                               grepl("Substantive", term) ~ "Substantive",
                               grepl("Personalization", term) ~ "Personalization",
                               grepl("SurrogationMP for another constituency", term) ~ "Territorial Surrogation",
                               grepl("PartisanSurrogationMP", term) ~ "Partisan Surrogation",
                               grepl("Responsiveness", term) ~ "Responsiveness",
                               grepl("Justification", term) ~ "Justification")) 

us_interaction_ests_group <- us_coefs_group %>% filter(grepl(":",term)) %>%
  mutate(group = case_when(grepl("GroupMid", term) ~ "Mid",
                           grepl("GroupHi", term) ~ "High"),
         dimension = case_when(grepl("Descriptive1", term) ~ "Descriptive (Partial)",
                               grepl("Descriptive2", term) ~ "Descriptive (Full)",
                               grepl("Substantive", term) ~ "Substantive",
                               grepl("Personalization", term) ~ "Personalization",
                               grepl("SurrogationRepresentative for another district", term) ~ "Territorial Surrogation",
                               grepl("PartisanSurrogation", term) ~ "Partisan Surrogation",
                               grepl("Responsiveness", term) ~ "Responsiveness",
                               grepl("Justification", term) ~ "Justification"))

de_interaction_ests_group <- de_coefs_group %>% filter(grepl(":",term)) %>%
  mutate(group = case_when(grepl("GroupMid", term) ~ "Mid",
                           grepl("GroupHi", term) ~ "High"),
         dimension = case_when(grepl("Descriptive1", term) ~ "Descriptive (Partial)",
                               grepl("Descriptive2", term) ~ "Descriptive (Full)",
                               grepl("Substantive", term) ~ "Substantive",
                               grepl("Personalization", term) ~ "Personalization",
                               grepl("SurrogationMP for another constituency in another", term) ~ "Territorial Surrogation", # Note this is only one of the DE Territorial surrogation dimensions
                               grepl("PartisanSurrogation", term) ~ "Partisan Surrogation",
                               grepl("Responsiveness", term) ~ "Responsiveness",
                               grepl("Justification", term) ~ "Justification")) %>%
  filter(!is.na(dimension))

uk_interaction_ests_group$country <- "UK"
us_interaction_ests_group$country <- "US"
de_interaction_ests_group$country <- "Germany"

interaction_group_out <- rbind(uk_interaction_ests_group, us_interaction_ests_group, de_interaction_ests_group) %>%
  mutate(dimension = factor(dimension, levels = rev(c("Substantive", "Descriptive (Partial)", "Descriptive (Full)", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")))) 

interaction_group_out %>%
  ggplot(aes(x = estimate, y = dimension, group = group, col = group)) + 
  geom_pointrange(aes(xmin = boot_low, xmax = boot_high), position = position_dodge(.2)) + 
  scale_color_manual("Factor group", values = c("black", "gray")) + 
  theme_bw() + 
  geom_vline(xintercept = 0, linetype = 2) + 
  xlab("Average Marginal Component Interaction Effect") + 
  ylab("") +
  facet_grid(~country)
ggsave(paste0("../output/plots/conjoint_amce_validation/amce_interaction_fa_group",attn,".pdf"), width = 12, height = 6)

sink(paste0("../output/tables/conjoint/amce_interaction_group",attn,".tex"))
interaction_group_out %>%
  mutate(est = paste0(round(estimate, 3), " (",round(std.error,3), ")")) %>%
  select(dimension, group, est, country) %>%
  pivot_wider(names_from = country,
              values_from = est) %>%
  arrange(desc(dimension), desc(group)) %>%
  rename(Dimension = dimension,
         `Factor Group`= group) %>%
  kableExtra::kable(format = "latex", booktabs = TRUE, escape = TRUE, 
                    align = "lllllll", linesep = "", 
                    caption = paste0("Average Marginal Component Interaction Effects (Binned Groups)."), 
                    label = paste0("amce_interaction_effects_binned",attn)) %>%
  kableExtra::kable_styling() %>%
  kableExtra::kable_styling(font_size = 10) %>%
  kableExtra::add_footnote(paste0("Sample sizes: Germany = ", nrow(de_long_clean),", UK = ",nrow(uk_long_clean),", US = ", nrow(us_long_clean),". The dependent variable is equal to 1 if an MP was selected as the better representative in a given comparison and 0 otherwise. Respondent-clustered standard errors in parentheses."), notation = "number") %>%
  print()
sink()

## Calculate subgroup AMCEs

calculate_conditional_amces <- function(coefs){

  interaction_terms <- grepl("\\:", names(coefs))
  c_amces <- rep(NA, length(interaction_terms))
  c_amces[!interaction_terms] <- coefs[!interaction_terms]
  
  constituent_terms <- unlist(lapply(strsplit(names(coefs)[interaction_terms], ":"), function(x) x[!grepl("GroupMid|GroupHi", x)]))
  
  for(i in 1:sum(interaction_terms)){
    
    c_amces[which(interaction_terms)[i]] <- coefs[which(interaction_terms)[i]] + coefs[names(coefs) == constituent_terms[i]]    
  }
  
  names(c_amces) <- names(coefs)
  c_amces
    
}

get_boot_ests <- function(model, boot_ests){
  
  camce_ests <- calculate_conditional_amces(coef(model))
  camce_boot <- apply(boot_ests,2, calculate_conditional_amces)
  camce_lo <- apply(camce_boot, 1, quantile, 0.025)
  camce_hi <- apply(camce_boot, 1, quantile, 0.975)
  camce <- data.frame(term = names(camce_ests), est = camce_ests, lo = camce_lo, hi = camce_hi)
  camce
    
}

camce_uk <- get_boot_ests(uk_conjoint_interaction, boot_uk_group)
camce_us <- get_boot_ests(us_conjoint_interaction, boot_us_group)
camce_de <- get_boot_ests(de_conjoint_interaction, boot_de_group)
camce_uk$country <- "UK"
camce_us$country <- "US"
camce_de$country <- "DE"

camce_combined <- tibble(rbind(camce_uk, camce_us, camce_de)) %>%
  mutate(group = case_when(grepl("GroupMid", term) ~ "Mid",
                           grepl("GroupHi", term) ~ "High",
                           !grepl("Group", term) ~ "Low"),
         group = factor(group, levels= c("Low", "Mid", "High")),
         dimension = case_when(grepl("Descriptive1", term) ~ "Descriptive\n (Partial)",
                               grepl("Descriptive2", term) ~ "Descriptive\n (Full)",
                               grepl("Substantive", term) ~ "Substantive",
                               grepl("Personalization", term) ~ "Personalization",
                               grepl("SurrogationRepresentative for another district|SurrogationMP for another constituency in another state", term) ~ "Territorial\n Surrogation",
                               grepl("PartisanSurrogation", term) ~ "Partisan\n Surrogation",
                               grepl("Responsiveness", term) ~ "Responsiveness",
                               grepl("Justification", term) ~ "Justification"),
         dimension = factor(dimension, levels = (c("Substantive", "Descriptive\n (Partial)", "Descriptive\n (Full)", "Territorial\n Surrogation", "Partisan\n Surrogation", "Justification", "Personalization", "Responsiveness")))) %>%
  filter(!is.na(dimension)) %>%
  filter(!(grepl("Group", term) & !grepl("\\:", term))) 

camce_combined %>%
  ggplot(aes(x = est, xmin = lo, xmax = hi, y = group)) + 
  geom_pointrange() + 
  theme_bw() + geom_vline(xintercept = 0, linetype = 2) + 
  scale_color_manual("Factor group", values = c("gray", "red", "black"))+ 
  facet_grid(dimension~country) + 
  ylab("Factor score group") + 
  xlab("Conditional Average Marginal Component Effect") + 
  theme(strip.background = element_blank(), strip.placement = "outside",
        legend.position = "bottom")

ggsave(paste0("../output/plots/conjoint_amce_validation/conditional_amce_fa_group",attn,".pdf"), width = 12, height = 9)


## F-tests for interactions

vars_to_sub <- c("Descriptive_FA_Group", 
                 "Substantive_FA_Group", 
                 "Personalization_FA_Group", 
                 "Surrogation_FA_Group", 
                 "Responsiveness_FA_Group",
                 "Justification_FA_Group")

i <- 1

out_list <- list()

for(i in 1:length(vars_to_sub)){

  tmp_form <- as.formula(paste0("mp_selected ~", gsub(paste0("\\* ",vars_to_sub[i]), "", as.character(f1_int)[3])))
  
  tmp_interaction_uk <- survey::svyglm(tmp_form, design = svydesign_out_uk)
  p_uk <- anova(tmp_interaction_uk, uk_conjoint_interaction)$p
  
  tmp_interaction_de <- survey::svyglm(tmp_form, design = svydesign_out_de)
  p_de <- anova(tmp_interaction_de, de_conjoint_interaction)$p
  
  tmp_interaction_us <- survey::svyglm(tmp_form, design = svydesign_out_us)
  p_us <- anova(tmp_interaction_us, us_conjoint_interaction)$p
  
  out <- data.frame(var = vars_to_sub[i], 
                    p_uk = p_uk, 
                    p_us = p_us,
                    p_de = p_de)
  
  out_list[[i]] <- out  

  }

f_tests <- do.call("rbind", out_list)
f_tests[,-1] <- round(f_tests[,-1],4)
f_tests$var <- gsub("_FA_Group","", f_tests$var)
f_tests[,-1] < 0.05
save(f_tests, file = paste0("../working/ftests",attn,".Rdata"))


## Interaction effects -- substantive representation x issue importance

f1_int_substantive_importance <- mp_selected ~ Descriptive + Substantive_Importance_Interaction + Personalization  + Surrogation + PartisanSurrogation + Responsiveness + Justification  

interaction_levels <- c("Incongruent_Low", "Incongruent_Mid", "Incongruent_High", 
                       "Congruent_Low", "Congruent_Mid", "Congruent_High")

uk_long_clean <- uk_long_clean %>% mutate(Substantive_Importance_High_Low = factor(case_when(Substantive_Importance %in% c(1:2) ~ "Low", 
                                                                                             Substantive_Importance %in% c(3:5) ~ "Mid",
                                                                                             Substantive_Importance %in% c(6:7) ~ "High"), 
                                                                                   levels = c("Low", "Mid", "High")),
                                          Substantive_Importance_Interaction = trimws(paste0(Substantive,"_", Substantive_Importance_High_Low)),
                                          Substantive_Importance_Interaction = ifelse(grepl("NA", Substantive_Importance_Interaction), NA, Substantive_Importance_Interaction),
                                          Substantive_Importance_Interaction = factor(Substantive_Importance_Interaction, levels = c(interaction_levels)))

us_long_clean <- us_long_clean %>% mutate(Substantive_Importance_High_Low = factor(case_when(Substantive_Importance %in% c(1:2) ~ "Low", 
                                                                                             Substantive_Importance %in% c(3:5) ~ "Mid",
                                                                                             Substantive_Importance %in% c(6:7) ~ "High"), 
                                                                                   levels = c("Low", "Mid", "High")),
                                          Substantive_Importance_Interaction = trimws(paste0(Substantive,"_", Substantive_Importance_High_Low)),
                                          Substantive_Importance_Interaction = ifelse(grepl("NA", Substantive_Importance_Interaction), NA, Substantive_Importance_Interaction),
                                          Substantive_Importance_Interaction = factor(Substantive_Importance_Interaction, levels = c(interaction_levels)))

de_long_clean <- de_long_clean %>% mutate(Substantive_Importance_High_Low = factor(case_when(Substantive_Importance %in% c(1:2) ~ "Low", 
                                                                                             Substantive_Importance %in% c(3:5) ~ "Mid",
                                                                                             Substantive_Importance %in% c(6:7) ~ "High"), 
                                                                                   levels = c("Low", "Mid", "High")),
                                          Substantive_Importance_Interaction = trimws(paste0(Substantive,"_", Substantive_Importance_High_Low)),
                                          Substantive_Importance_Interaction = ifelse(grepl("NA", Substantive_Importance_Interaction), NA, Substantive_Importance_Interaction),
                                          Substantive_Importance_Interaction = factor(Substantive_Importance_Interaction, levels = c(interaction_levels)))


amces_uk <- cj(uk_long_clean, f1_int_substantive_importance, id = ~ResponseId, weights = uk_long_clean$wgt)
amces_us <- cj(us_long_clean, f1_int_substantive_importance, id = ~ResponseId, weights = us_long_clean$wgt)
amces_de <- cj(de_long_clean, f1_int_substantive_importance, id = ~ResponseId, weights = de_long_clean$wgt)

amces_us <- amces_us %>% 
  mutate(feature = gsub("PartisanSurrogation","Partisan Surrogation",feature),
         feature = gsub("^Surrogation", "Territorial Surrogation", feature),
         feature = gsub("Substantive_Importance_Interaction", "Substantive", feature),
         level = as.character(level),
         level = case_when(level == "Congruent_High" ~ "Congruent (High Salience)",
                           level == "Congruent_Low" ~ "Congruent (Low Salience)",
                           level == "Congruent_Mid" ~ "Congruent (Mid Salience)",
                           level == "Incongruent_High" ~ "Incongruent (High Salience)",
                           level == "Incongruent_Low" ~ "Incongruent (Low Salience)",
                           level == "Incongruent_Mid" ~ "Incongruent (Mid Salience)",
                           TRUE ~ level),
         level = factor(level, levels = level),
         feature = factor(feature, levels = c("Substantive", "Descriptive", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")))

amces_uk <- amces_uk %>% 
  mutate(feature = gsub("PartisanSurrogation","Partisan Surrogation",feature),
         feature = gsub("^Surrogation", "Territorial Surrogation", feature),
         feature = gsub("Substantive_Importance_Interaction", "Substantive", feature),
         level = as.character(level),
         level = case_when(level == "Congruent_High" ~ "Congruent (High Salience)",
                           level == "Congruent_Low" ~ "Congruent (Low Salience)",
                           level == "Congruent_Mid" ~ "Congruent (Mid Salience)",
                           level == "Incongruent_High" ~ "Incongruent (High Salience)",
                           level == "Incongruent_Low" ~ "Incongruent (Low Salience)",
                           level == "Incongruent_Mid" ~ "Incongruent (Mid Salience)",
                           TRUE ~ level),
         level = factor(level, levels = level),
         feature = factor(feature, levels = c("Substantive", "Descriptive", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")))

amces_de <- amces_de %>% 
  mutate(feature = gsub("PartisanSurrogation","Partisan Surrogation",feature),
         feature = gsub("^Surrogation", "Territorial Surrogation", feature),
         feature = gsub("Substantive_Importance_Interaction", "Substantive", feature),
         level = as.character(level),
         level = case_when(level == "Congruent_High" ~ "Congruent (High Salience)",
                           level == "Congruent_Low" ~ "Congruent (Low Salience)",
                           level == "Congruent_Mid" ~ "Congruent (Mid Salience)",
                           level == "Incongruent_High" ~ "Incongruent (High Salience)",
                           level == "Incongruent_Low" ~ "Incongruent (Low Salience)",
                           level == "Incongruent_Mid" ~ "Incongruent (Mid Salience)",
                           TRUE ~ level),
         level = factor(level, levels = level),
         feature = factor(feature, levels = c("Substantive", "Descriptive", "Territorial Surrogation", "Partisan Surrogation", "Justification", "Personalization", "Responsiveness")))

xlims <- range(c(amces_uk$lower,amces_uk$upper, amces_us$lower, amces_us$upper, amces_de$lower, amces_de$upper), na.rm = TRUE)


amces_de$country <- "Germany"
amces_uk$country <- "UK"
amces_us$country <- "US"

amces_all <- rbind(amces_de, amces_uk, amces_us)

levels(amces_all$level) <- gsub("district|constituency","district", levels(amces_all$level))
levels(amces_all$level) <- gsub("MP","Representative", levels(amces_all$level))
levels(amces_all$level) <- gsub("Representative","Rep.", levels(amces_all$level))
amces_all$feature <- gsub("Territorial ", "Territorial \n", amces_all$feature)
amces_all$feature <- gsub("Partisan ", "Partisan \n", amces_all$feature)

amces_all$feature <- factor(amces_all$feature, levels = c("Substantive", "Descriptive", "Territorial \nSurrogation", "Partisan \nSurrogation", "Justification", "Personalization", "Responsiveness"))

ggplot(amces_all, aes(x = estimate, xmax = upper, xmin = lower, y = level, group = feature)) + 
  geom_point() + geom_pointrange() + 
  facet_grid(feature ~ country, scales = "free_y" ) + 
  theme_bw() + geom_vline(xintercept = 0, linetype = 2) + 
  ylab("") + xlab("AMCE") + 
  theme(strip.background = element_blank(), strip.placement = "outside")

ggsave(paste0("../output/plots/conjoint_amce/conjoint_amce_combined_interaction_salience",attn,".pdf"), width = 16, height = 7.5)


salience_tab <- amces_all %>%
  mutate(est = paste0(round(estimate, 3), " (", round(std.error,3),")"),
         est = ifelse(grepl("NA", est), "-", est)) %>%
  filter(est != "-") %>%
  select(feature, level, est, country) %>%
  pivot_wider(values_from = est,
              names_from = country) %>%
  mutate(feature = gsub("\\n", "", as.character(feature)),
         level = trimws(as.character(level))) %>%
  rename(Dimension = feature,
         Level = level) %>%
  kableExtra::kable(format = "latex", booktabs = TRUE, escape = TRUE, 
                    align = "lllllll", linesep = "", 
                    caption = paste0("Average Marginal Component Effects by Issue Salience."), 
                    label = paste0("amce_by_issue_salience")) %>%
  kableExtra::kable_styling() %>%
  kableExtra::kable_styling(font_size = 10) %>%
  kableExtra::add_footnote(paste0("Sample sizes: Germany = ", nrow(de_long_clean),", UK = ",nrow(uk_long_clean),", US = ",nrow(us_long_clean),". The dependent variable is equal to 1 if an MP was selected as the better representative in a given comparison and 0 otherwise. Respondent-clustered standard errors in parentheses."), notation = "number")

sink(paste0("../output/tables/conjoint/conjoint_amce_combined_issue_salience_conditional",attn,".tex"))
print(salience_tab)
sink()

# Calculate effect sizes for salient and non-salient issues
range(amces_all$estimate[amces_all$level == "Congruent (High Salience)"] - amces_all$estimate[amces_all$level == "Incongruent (High Salience)"])

range(amces_all$estimate[amces_all$level == "Congruent (Low Salience)"] - amces_all$estimate[amces_all$level == "Incongruent (Low Salience)"])

